library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.2 --
## v ggplot2 3.4.3     v purrr   0.3.4
## v tibble  3.2.1     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.5.1
## v readr   2.1.3     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(metaumbrella)
collapsunique <- function(x) paste(unique(sort(x)), collapse = ", ")

source("1b - cct_level_cam_scoring_140424.R")
## New names:
## New names:
## Joining, by = c("first_author_meta", "year_meta")
## Joining, by = c("first_author_meta", "year_meta", "AMS_2_protocol",
## "AMS_4_search", "AMS_5_screening_duplicate", "AMS_6_extraction_duplicate",
## "AMS_9_RoB")
## Joining, by = c("Reference", "ref_year", "doi", "first_author_meta",
## "year_meta", "status", "amstar_done", "acronym", "id_row_meta", "therap",
## "AMS_1_PICO", "AMS_2_protocol", "AMS_3_selection", "AMS_4_search",
## "AMS_5_screening_duplicate", "AMS_6_extraction_duplicate", "AMS_7_excluded",
## "AMS_8_description", "AMS_9_RoB", "AMS_10_funding_primary", "AMS_11_stat",
## "AMS_12_MA_RoB", "AMS_13_interpreting_RoB", "AMS_14_heterogeneity",
## "AMS_15_publication_bias", "AMS_16_CoI", "Link_MA", "...28", "...29")
## Joining, by = c("first_author_meta", "year_meta", "intervention_general",
## "author", "year", "trial_id_supp")
## Joining, by = c("first_author_meta", "year_meta", "intervention_general",
## "author", "year", "trial_id_supp", "paper", "id_trial_meta", "id_trial")
## * `` -> `...40`
## * `` -> `...41`
## * `` -> `...42`
## * `` -> `...43`
## * `` -> `...44`
## * `` -> `...45`
source("2 - cct_level_homogeneization_140424.R")
## New names:
## * `...28` -> `...105`
## * `...29` -> `...106`
source("3 - ma_level_scoring140424.R")
## 
## Attaching package: 'Hmisc'
## 
## The following objects are masked from 'package:dplyr':
## 
##     src, summarize
## 
## The following objects are masked from 'package:base':
## 
##     format.pval, units
## 
## `summarise()` has grouped output by 'paper', 'intervention_general', 'intervention_precise'. You can override using the `.groups` argument.
chemin = "C:/Users/Corentin Gosling/drive_gmail/Recherche/Article 2 - Base de Donnees/7 - Data analysis/data/"
# chemin = "H:/Mon Drive/Recherche/Article 2 - Base de Donnees/7 - Data analysis/data/"

dcct = readxl::read_excel(paste0(chemin,
                                "cct_level_homogeneized_scored.xlsx")) 
dma = readxl::read_excel(paste0(chemin,
                                "ma_level_homogeneized_scored.xlsx"))
dcct$measure[dcct$measure == "SMC"] = "SMD"
dcct$measure[dcct$measure == "MC"] <- "MD"

dcct = dcct %>% filter(intervention_type == "Complementary")
dcct$factor = with(dcct, paste0(paper, "_", intervention_general, "_", outcome_general))
dma$Factor = with(dma, paste0(paper, "_", intervention_general, "_", outcome_general))
# dma[dma$paper == "Fraguas (2019)", c("Factor", "mean_age")]
# View(dma %>% filter(Factor == "Martins (2021)_OXYT_Restricted/repetitive behaviors") %>% select(Factor, mean_age, mean_age_q10, mean_age_q15, mean_age_q20, mean_age_q25))



dat_error = view.errors.umbrella(dcct)
## ERROR:
## - Some repeated studies (author and year) in the same factor do not have any 'multiple_es' value.
## - Study with several multiple_es values. Please, specify only one of either 'groups' or 'outcomes'. 
## WARNING:
## - A numeric cell in the column 'value' contains a ',' that was converted to a '.' (this is only a warning, not an error).
dcct[dat_error[dat_error$column_type_errors == "ERROR", "row_index"], "multiple_es"] <- "outcomes"
dat_error = view.errors.umbrella(dcct)
## WARNING:
## - A numeric cell in the column 'value' contains a ',' that was converted to a '.' (this is only a warning, not an error).
## - A numeric cell in the column 'ci_lo' contains a ',' that was converted to a '.' (this is only a warning, not an error).
## - A numeric cell in the column 'ci_up' contains a ',' that was converted to a '.' (this is only a warning, not an error).
## - A numeric cell in the column 'mean_cases' contains a ',' that was converted to a '.' (this is only a warning, not an error).
## - A numeric cell in the column 'sd_cases' contains a ',' that was converted to a '.' (this is only a warning, not an error).
## - A numeric cell in the column 'mean_controls' contains a ',' that was converted to a '.' (this is only a warning, not an error).
## - A numeric cell in the column 'sd_controls' contains a ',' that was converted to a '.' (this is only a warning, not an error).
## - The confidence interval is not symmetric around the effect size (this is only a warning, not an error).
## - Some but not all effect sizes are reversed for a factor. Check this is what you want (this is only a warning, not an error).
# rio::export(dcct, paste0("C:/Users/Corentin Gosling/drive_gmail/Recherche/Article 2 - Base de Donnees/7 - Data analysis/data/", "dat-CAM-analysis.txt"))

SX. Primary analysis

SX-1. Data analysis

# library(tidyverse)
# library(metaumbrella)
# # dcct = read.delim(paste0("C:/Users/Corentin Gosling/drive_gmail/Recherche/Article 2 - Base de Donnees/7 - Data analysis/data/", "dat-CAM-analysis.txt"))
# res = umbrella(dcct %>%
#                  filter(first_author_meta %in% c("Zhou", "Salazar de Pablo")),

res = umbrella(dcct,
               verbose=FALSE,
               mult.level = TRUE, method.var = "REML", 
               r = 0.8, pre_post_cor = 0.5)
res_m_cred = summary(res)
res_m_cred$Outcome = t(do.call(cbind, stringr::str_split(res_m_cred$Factor, "_")))[, 3]
res_m_cred$Intervention = t(do.call(cbind, stringr::str_split(res_m_cred$Factor, "_")))[, 2]
res_m_cred$"Meta-review" = t(do.call(cbind, stringr::str_split(res_m_cred$Factor, "_")))[, 1]
res_m = left_join(res_m_cred, dma, by="Factor") 

# View(res_m %>%
#        select(Factor, starts_with("available"), IN_meta, amstar.x, amstar.y, active_controls, mean_age))

SX-1. Certainty of evidence

Scoring components

RoB

res_m$down_rob = with(res_m,
    ifelse(rob > 75, 0, 
           ifelse(rob > 50, 1, 2)))
# View(res_m[, c("down_rob", "rob", "Factor")])

Inconsistency

# res_m$min_es = res_m$max_es = res_m$n_agreement = NA
# i = 0
# res_m$value = as.numeric(as.character(res_m$value))
# for (fact in unique(res_m$Factor)) {
#   i = i+1
#   # print(fact)
#   # print(i)
#   if (fact != res_m$Factor[i]) {
#     stop(paste0(fact, " != ", res_m$Factor[i], "(", i, ")"))
#   }
#   dat_sub = NA
#   dat_sub = res[[fact]]$x
#   values=pooled_es=NA
#   if (res_m$measure[i] == "G") {
#     values = dat_sub$value
#     pooled_es = res_m[res_m$Factor == fact, ]$value
#   } else {
#     values = log(dat_sub$value)
#     pooled_es = log(res_m[res_m$Factor == fact, ]$value)
#   }
#   res_m$min_es[i] = min(values)
#   res_m$max_es[i] = max(values)
#   res_m$n_agreement[i] = sum(sign(values) == sign(pooled_es)) / length(dat_sub$value)
# }
# res_m$down_het = with(res_m,
#     ifelse((I2 >= 50 & n_agreement < 0.90), 2,
#            ifelse(I2 >= 25 & n_agreement < 0.90, 1, 0)))

# View(res_m[, c("n_studies", "min_es", "max_es", "n_agreement", "I2", "down_het", "Factor")])

eq_range_or = c(0.8, 1.25)
eq_range_g = c(-0.1, 0.1)
meas_G = res_m$measure == "G"
meas_OR = res_m$measure %in% c("OR", "RR")

low_ci_pos = (meas_G & CI_lo_value > 0) | (meas_OR & CI_lo_value > 1)
low_ci_neg_range = (meas_G & CI_lo_value < 0) & (meas_G & CI_lo_value > eq_range_g[1]) |
                   (meas_OR & CI_lo_value < 1) & (meas_OR & CI_lo_value > eq_range_or[1])
low_ci_out_range = (meas_G & CI_lo_value < eq_range_g[1]) |
                   (meas_OR & CI_lo_value < eq_range_or[1])

low_pi_pos = (meas_G & PI_lo_value > 0) | (meas_OR & PI_lo_value > 1)
low_pi_neg_range = (meas_G & PI_lo_value < 0 & PI_lo_value > eq_range_g[1]) |
                   (meas_OR & PI_lo_value < 1) & (meas_OR & PI_lo_value > eq_range_or[1])
low_pi_out_range = (meas_G & PI_lo_value < eq_range_g[1]) |
                   (meas_OR & PI_lo_value < eq_range_or[1]) 

up_ci_neg = (meas_G & CI_up_value < 0) | (meas_OR & CI_up_value < 1)
up_ci_pos_range = (meas_G & CI_up_value > 0 & CI_up_value < eq_range_g[2]) |
                  (meas_OR & CI_up_value > 1) & (meas_OR & CI_up_value < eq_range_or[2])
up_ci_out_range = (meas_G & CI_up_value > eq_range_g[2]) |
                  (meas_OR & CI_up_value > eq_range_or[2])  

up_pi_neg = (meas_G & PI_up_value < 0) | (meas_OR & PI_up_value < 1)
up_pi_pos_range = (meas_G & PI_up_value > 0) & (meas_G & PI_up_value < eq_range_g[2]) |
                  (meas_OR & PI_up_value > 1) & (meas_OR & PI_up_value < eq_range_or[2])
up_pi_out_range = (meas_G & PI_up_value > eq_range_g[2]) |
                  (meas_OR & PI_up_value > eq_range_or[2]) 

# res_m$DOWN_HET = ifelse((low_ci_pos & low_pi_pos) | (up_ci_neg & up_pi_neg), 0,
#   ifelse((low_ci_neg_range & low_pi_neg_range) | (up_ci_pos_range & up_pi_pos_range), 0,
#   ifelse((low_ci_out_range & low_pi_out_range) | (up_ci_out_range & up_ci_out_range), 0,
#   ifelse((low_ci_pos & low_pi_neg_range) | (up_ci_neg & up_pi_pos_range), 1,
#   ifelse((low_ci_neg_range & low_pi_out_range) | (up_ci_pos_range & up_pi_out_range), 1,
#   ifelse((low_ci_pos & low_pi_out_range) | (up_ci_neg & up_ci_out_range), 2,
#   ifelse((low_ci_neg_range & low_pi_out_range) & (up_ci_pos_range & up_pi_out_range), 2, 999)))))))
res_m$down_het = 
  ifelse((low_ci_neg_range & low_pi_out_range) & (up_ci_pos_range & up_pi_out_range), 2, 
  ifelse((low_ci_pos & low_pi_out_range) | (up_ci_neg & up_ci_out_range), 2,
  ifelse((low_ci_neg_range & low_pi_out_range) | (up_ci_pos_range & up_pi_out_range), 1, 
  ifelse((low_ci_pos & low_pi_neg_range) | (up_ci_neg & up_pi_pos_range), 1,
  ifelse((low_ci_pos & low_pi_pos) | (up_ci_neg & up_pi_neg), 0,
  ifelse((low_ci_neg_range & low_pi_neg_range) | (up_ci_pos_range & up_pi_pos_range), 0,
  ifelse((low_ci_out_range & low_pi_out_range) | (up_ci_out_range & up_ci_out_range), 0,
  999)))))))

rows_smd = which(res_m$measure == "G")
rows_rr = which(res_m$measure != "G")

row_imp_pi = which(is.na(res_m$PI_lo_value) | is.na(res_m$PI_up_value))
res_m$down_het[row_imp_pi] <- 2


# View(cbind(res_m[, c("PI_lo_value", "CI_lo_value", "CI_up_value", "PI_up_value", 
#                "measure", "DOWN_HET", "n_studies")],
#                 low_ci_pos, low_ci_neg_range, low_ci_out_range,
#                low_pi_pos, low_pi_neg_range, low_pi_out_range,
#                up_ci_neg, up_ci_pos_range, up_ci_out_range,
#                up_pi_neg, up_pi_pos_range, up_pi_out_range)[rows_rr,])

Indirectness

res_m$down_ind =  with(res_m,
    ifelse(grepl("Mixed", res_m$PICO_amstar_ID, fixed = TRUE), 1,
    ifelse(res_m$available_control < 0.75 , 1, 0)
      # (active_controls >= 80 & active_controls <= 20) |
    ))

Imprecision

res_m$cross_low_high = (meas_G & (res_m$CI_lo_value < 0 & res_m$CI_up_value > 0.8) |
                          (res_m$CI_lo_value < -0.8 & res_m$CI_up_value > 0)) | 
  (meas_OR & (res_m$CI_lo_value < 1 & res_m$CI_up_value > 5) |
                          (res_m$CI_lo_value < 0.2 & res_m$CI_up_value > 1))

res_m$n_fail_detect_small_effects = 
  res_m$n_cases < 394 | res_m$n_controls < 394
res_m$n_fail_detect_large_effects = 
  res_m$n_cases < 64 | res_m$n_controls < 64

res_m$down_imp = ifelse(
  (res_m$cross_low_high & res_m$n_fail_detect_small_effects) |
    res_m$n_fail_detect_large_effects, 2,
  ifelse(res_m$cross_low_high | res_m$n_fail_detect_small_effects, 1, 0))

# View(res_m[, c( "down_imp", "ci_lo", "ci_up", "cross_low_high", "n_cases", "n_controls", "Factor")])

Publication bias

res_m$down_pubbias = ifelse(
  as.numeric(res_m$egger_p) < 0.10 |
  res_m$n_studies < 5 | 
  is.na(as.numeric(res_m$egger_p)), 1, 0)
## Warning in ifelse(as.numeric(res_m$egger_p) < 0.1 | res_m$n_studies < 5 | : NAs
## introduced by coercion

## Warning in ifelse(as.numeric(res_m$egger_p) < 0.1 | res_m$n_studies < 5 | : NAs
## introduced by coercion

Overal scoring

res_m$down_rob[is.na(res_m$down_rob)] <- 2
res_m$down_het[is.na(res_m$down_het)] <- 2
res_m$down_ind[is.na(res_m$down_ind)] <- 2
res_m$down_imp[is.na(res_m$down_imp)] <- 2
res_m$down_pubbias[is.na(res_m$down_pubbias)] <- 2

res_m$GRADE_num = rowSums(res_m[, 
        c("down_rob", "down_het", "down_ind", "down_imp", "down_pubbias")]) 
  
res_m$GRADE = with(res_m, ifelse(
  GRADE_num == 0, "High", 
  ifelse(GRADE_num %in% 1, "Moderate", 
         ifelse(GRADE_num %in% 2, "Low", 
                ifelse(GRADE_num >= 3, "Very low", NA)))))

res_m$Age = factor(res_m$Age, 
                   levels=c('< 6 yo', '6-12 yo', '13-19 yo', '>= 20 yo'))
res_m$GRADE = factor(res_m$GRADE, 
                   levels=c('High', 'Moderate', 'Low', 'Very low'))

res_m = res_m %>%
  arrange(Age, Intervention, GRADE) 

View(res_m %>%
         filter(GRADE == "Moderate") %>% 
       select(PICO_amstar_ID, IN_meta, n_studies, eG, eG_CI, 
              down_rob, down_het, down_ind, down_imp, down_pubbias, 
              I2, n_cases, n_controls, PI_eG))

colSums(res_m[, 
        c("down_rob", "down_het", "down_ind", 
          "down_imp", "down_pubbias")] == 2)
##     down_rob     down_het     down_ind     down_imp down_pubbias 
##           92           94            0          135            0
colSums(res_m[, 
        c("down_rob", "down_het", "down_ind", 
          "down_imp", "down_pubbias")] == 1)
##     down_rob     down_het     down_ind     down_imp down_pubbias 
##           60           36           83          112          169
res_m$paper = gsub("\\(child\\) ", "", res_m$paper)
res_m$paper = gsub("\\(adult\\) ", "", res_m$paper)
res_m$n_paper = length(unique(paste0(res_m$paper)))
interventions = sort(unique(res_m$intervention_general))
  
Mind_body_medicine = c("MUSIC", "SENS", "PHYS", "rTMS", "tCDS")

therapies_energetiques = c("ACUP");

Natural_Product_Based_Therapies  = c(
  "DIET", "FOLI", "HERB", "L-CARNIT", "L-CARNO", 
  "MELAT", "NAC", "OXYT", "PROB", "PUFA", 
  "SECRET", "SULFO", "V1a-RA", "VIT-D")
systemes_medicaux_alternatifs = c("AAI")

outcomes = sort(unique(res_m$outcome_general)); interventions; outcomes
##  [1] "AAI"      "ACUP"     "DIET"     "FOLI"     "HERB"     "L-CARNIT"
##  [7] "L-CARNO"  "MELAT"    "MUSIC"    "NAC"      "OXYT"     "PHYS"    
## [13] "PROB"     "PUFA"     "rTMS"     "SECRET"   "SENS"     "SULFO"   
## [19] "tDCS"     "V1a-RA"   "VIT-D"
##  [1] "Acceptability"                   "Adaptative Functioning"         
##  [3] "ADHD symptoms"                   "Adverse events"                 
##  [5] "Anxiety"                         "Disruptive behaviors"           
##  [7] "Global cognition (IQ)"           "Language (Expressive skills)"   
##  [9] "Language (Overall skills)"       "Language (Receptive skills)"    
## [11] "Mood related symptoms"           "Overall ASD symptoms"           
## [13] "Quality of life"                 "Restricted/repetitive behaviors"
## [15] "Sensory Profile"                 "Sleep quality"                  
## [17] "Sleep quantity"                  "Social-communication"           
## [19] "Tolerability"

Primary analysis

SX-1. Synthesis - table

res_p = res_m %>% filter(IN_meta == 1)

DT::datatable(res_p)
# res_p %>%
#         filter(GRADE %in% c("Low", "Moderate") ) %>%
#         select(PICO_amstar, GRADE, Age, eG, p_value) %>%
#         arrange(Age, GRADE)
# res_m %>%
#         filter(Intervention == "OXYT" & Outcome == "Restricted/repetitive behaviors") %>%
#         select(PICO_amstar, GRADE, Age, eG, p_value, paper) %>%
#         arrange(Age, GRADE)

SX-1. Synthesis - Plots

Low or Moderate GRADE rating

dat_low_or_higher = res_p %>% filter(GRADE %in% c("Low", "Moderate")) %>%
  mutate(age_pre = case_when(Age == "< 6 yo" ~ "Pre-school (<6 years old)",
                             Age == "6-12 yo" ~ "School-age (6-12 years old)",
                             Age == "13-19 yo" ~ "Adolescents (13-19 years old)",
                             Age == ">= 20 yo" ~ "Adults (>19 years old)"))
metaumbrella::forest(dat_low_or_higher,
       layout = "RevMan5",
       subgroup = "age_pre",
       subgroup.name = "",
       leftcols = c("Intervention", "Outcome", "GRADE", 
                    "n_studies",  "rob_num",  "I2", "effect.ci"),
       leftlabs = c("Intervention", "Outcome", "GRADE",
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Low or Moderate GRADE"
)

Core symptoms

res_core = res_p %>% filter(Outcome %in% ASD_symptoms)
metaumbrella::forest(res_core,
       layout = "RevMan5",
       subgroup = "Outcome",
       subgroup.name = "",
       leftcols = c("Intervention", "Age", "GRADE", 
                    "n_studies",  "rob_num",  "I2", "effect.ci"),
       leftlabs = c("Intervention", "Age", "GRADE",
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Core ASD symptoms"
)

Safety

forest(res_p %>% filter(Outcome %in% safety),
       title = "equivalent Risk Ratio (eRR)",
       xlim = c(0.05, 20),
       layout = "RevMan5",
       subgroup = "Outcome",
       xlab = "Equivalent Risk Ratio (eRR)",
       subgroup.name = "",
       leftcols = c("Intervention", "Age", "GRADE", 
                    "n_studies",  "rob_num",  "I2", "effect.ci"),
       leftlabs = c("Intervention", "Age", "GRADE", 
                    "n-studies", "Low\nRoB",   "I²", "eRR + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Safety"
)

Functioning

forest(res_p %>% filter(Outcome %in% functioning),
       layout = "RevMan5",
       subgroup = "Outcome",
       subgroup.name = "",
       leftcols = c("Intervention", "Age", "GRADE", 
                    "n_studies",  "rob_num",  "I2", "effect.ci"),
       leftlabs = c("Intervention", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Daily functioning"
)

Comorbidities

forest(res_p %>% filter(Outcome %in% comorbidities),
       layout = "RevMan5",
       subgroup = "Outcome",
       subgroup.name = "",
       leftcols = c("Intervention", "Age", "GRADE", 
                    "n_studies",  "rob_num",  "I2", "effect.ci"),
       leftlabs = c("Intervention", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Psychiatric comorbidities"
)

SX2. Overlapping

Overall picture

res_m$rob_num = paste0(round(res_m$rob), "%")

set.seed(123)
dat <- data.frame(
  intervention = rep(letters[1:10], each = 2),
  lower_bound = c(-abs(round(rnorm(20, 10, 3)))),
  upper_bound = c(abs(round(rnorm(20, 10, 3))))
)

# Function to calculate interval overlap
interval_overlap_percentage <- function(lower1, upper1, lower2, upper2) {
  overlap = max(0, min(upper1, upper2) - max(lower1, lower2))
  total_length = (upper1 - lower1) + (upper2 - lower2) - overlap
  percentage_overlap = (overlap / total_length) * 100
  return(percentage_overlap)
}
res_over <- res_m %>%
  group_by(PICO_amstar) %>%
  filter(n() > 1) %>%
  ungroup() %>%
  mutate(n_SRMA = length(unique(paper))) %>%
  group_by(PICO_amstar) %>%
  summarise(
    n = n(),
    n_SRMA = unique(n_SRMA),
    paper = collapsunique(paper),
    min_es = min(eG),
    max_es = max(eG),
    dif_es = max(eG) - min(eG),
    max_grade = head(sort(GRADE), 1),
    min_grade = tail(sort(GRADE), 1),
    prop_sig = sum(as.numeric(p_value) < 0.05) / n(),
    avg_percentage_overlap = mean(map_dbl(combn(n(), 2, simplify = FALSE), function(idx) {
    interval_overlap_percentage(
      ci_lo_g[idx[1]], ci_up_g[idx[1]], 
      ci_lo_g[idx[2]], ci_up_g[idx[2]])
  })))

propC = function(x, R=0) {
  x_non_na = x[!is.na(x)]
  round(sum(x_non_na)/length(x) * 100, R)
}

DT::datatable(res_over)
# identify, for each 

Mind-Body Medicine

MUSIC

forest(res_m %>% filter(Intervention %in% "MUSIC"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR (key items)", "AMSTAR (Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (MUSIC)")

SENS

forest(res_m %>% filter(Intervention %in% "SENS"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR (key items)", "AMSTAR (Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (SENS)")

PHYS

forest(res_m %>% filter(Intervention %in% "PHYS"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR (key items)", "AMSTAR (Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (PHYS)")

rTMS

forest(res_m %>% filter(Intervention %in% "rTMS"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR (key items)", "AMSTAR (Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (rTMS)")

tCDS

forest(res_m %>% filter(Intervention %in% "tDCS"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (tDCS)")

Energy Medicine

ACUP

forest(res_m %>% filter(Intervention %in% "ACUP"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (ACUP)")

Natural Product

Hormones

MELAT
forest(res_m %>% filter(Intervention %in% "MELAT"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (MELAT)")

V1a-RA
forest(res_m %>% filter(Intervention %in% "V1a-RA"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (V1a-RA)")

SECRET
forest(res_m %>% filter(Intervention %in% "SECRET"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (SECRET)")

OXYT
forest(res_m %>% filter(Intervention %in% "SECRET"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (SECRET)")

Supplementation

Natural_Product_Based_Therapies = c( ““,”“,”“,”“,”“,”“,”“,”“,”“,”“,”“,”“,”“,”“)

FOLI
forest(res_m %>% filter(Intervention %in% "FOLI"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (FOLI)")

HERB
forest(res_m %>% filter(Intervention %in% "HERB"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (HERB)")

L-CARNIT
forest(res_m %>% filter(Intervention %in% "L-CARNIT"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (L-CARNIT)")

L-CARNO
forest(res_m %>% filter(Intervention %in% "L-CARNO"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (L-CARNO)")

NAC
forest(res_m %>% filter(Intervention %in% "NAC"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (NAC)")

PROB
forest(res_m %>% filter(Intervention %in% "PROB"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (L-CARNIT)")

SULFO
forest(res_m %>% filter(Intervention %in% "SULFO"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (SULFO)")

PUFA
forest(res_m %>% filter(Intervention %in% "PUFA"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (PUFA)")

VIT-D
forest(res_m %>% filter(Intervention %in% "VIT-D"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (VIT-D)")

Diets

forest(res_m %>% filter(Intervention %in% "DIET"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (DIET)")

Whole Medical Systems

AAI
forest(res_m %>% filter(Intervention %in% "AAI"),
       layout = "RevMan5", subgroup = "Outcome",subgroup.name = "",
       leftcols = c("paper", "amstar.x", "amstar_rank", "Age", "GRADE", 
                    "n_studies",  "rob",  "I2", "effect.ci"),
       leftlabs = c("paper", "AMSTAR\n(key items)", "AMSTAR\n(Rank)", "Age", "GRADE", 
                    "n-studies", "Low RoB",   "I²", "eSMD + 95% CI"),
       smlab = "Overlapping (AAI)")

SX. Writing

Characteristics of the meta-analytic reports

paste0("The ", unique(res_m$n_paper), " retained reports explored the effects of ", 
       length(unique(res_m$intervention_general)), " intervention types, on at least one of the ", 
       length(unique(res_m$outcome_general)), " possible outcomes.",
       "These reports generated a total of ", nrow(res_m), " meta-analyses",
       " that synthesized the evidence from more than ",
       length(unique(paste0(dcct$author, dcct$year))), " CCTs, and that explored ", length(unique(res_m$PICO_amstar)), " unique PICOs.")
## [1] "The 53 retained reports explored the effects of 21 intervention types, on at least one of the 19 possible outcomes.These reports generated a total of 248 meta-analyses that synthesized the evidence from more than 206 CCTs, and that explored 147 unique PICOs."
res_year = res_m %>%
  group_by(paper) %>%
  slice(1)
paste0("These reports were, in their vast majority (", 
       propC(res_year$year_meta_num >= 2018), "%) published after 2018.")
## [1] "These reports were, in their vast majority (87%) published after 2018."
dat_pico_overlap = res_m %>%
  group_by(PICO_amstar, Intervention, Outcome, Age) %>%
  summarise(n=n()) %>%
  filter(n > 1)
## `summarise()` has grouped output by 'PICO_amstar', 'Intervention', 'Outcome'.
## You can override using the `.groups` argument.
paste0("A total of ", nrow(dat_pico_overlap), " PICOs (", round(nrow(dat_pico_overlap)/length(unique(res_m$PICO_amstar)) * 100), "%) was explored in several overlapping meta-analyses. The PICOs with the highest number of overlapping meta-analyses were investigating the effects of ", 
        unique(dat_pico_overlap[dat_pico_overlap$n >= 6, ]$Intervention), " on ", tolower(paste(paste0(paste0(unique(dat_pico_overlap[dat_pico_overlap$n >= 6, ]$Outcome)), " (n=",
paste0(unique(dat_pico_overlap[dat_pico_overlap$n >= 6, ]$n)), ")")
, collapse=", ")))
## [1] "A total of 53 PICOs (36%) was explored in several overlapping meta-analyses. The PICOs with the highest number of overlapping meta-analyses were investigating the effects of PUFA on adhd symptoms (n=7), disruptive behaviors (n=8), restricted/repetitive behaviors (n=7), social-communication (n=8)"
# paste0("The mean number of overlapping meta-analyses on the same PICO was ",
#       round(mean(dat_pico$n), 2), " (minimum = ", 
#       min(dat_pico$n), ", maximum = ",
#       max(dat_pico$n), ").")

# paste0(length(unique(res_m$intervention_general)), " intervention types (", sort(collapsunique(res_m$intervention_general)), ")")
# paste0(length(unique(res_m$outcome_general)), " outcomes (", sort(collapsunique(res_m$outcome_general)), ")")

Primary

res_interventions = res_p %>%
  group_by(Intervention) %>%
  select(Outcome) %>%
  distinct() %>%
  summarise(n = n(),
            outcome = collapsunique(Outcome)) %>%
  arrange(desc(n))
## Adding missing grouping variables: `Intervention`
nrowInter = nrow(res_interventions)

table(res_p$amstar_rank)
## 
## Critically low           High            Low 
##             34             24             90
dat_amstar = as.data.frame(table(res_p %>% group_by(paper) %>% slice(1) %>% ungroup%>% select(amstar_rank)))

paste0("Thus, ", length(unique(res_m$PICO_amstar)), " meta-analyses, derived from ", length(unique(res_p$paper)), " reports, were finally included in the primary analysis after discarding overlapping meta-analyses.")
## [1] "Thus, 147 meta-analyses, derived from 31 reports, were finally included in the primary analysis after discarding overlapping meta-analyses."
paste0("According to the AMSTAR-2 scoring, ",
        dat_amstar[dat_amstar$Var1=="High", "Freq"],
        " of these ", length(unique(res_p$paper)), 
        " reports were of high quality, ",
        dat_amstar[dat_amstar$Var1=="Low", "Freq"], 
        " of low quality, and ",
        dat_amstar[dat_amstar$Var1=="Critically low", "Freq"],
        " of critically low quality. ",
        dat_amstar[dat_amstar$Var1=="Moderate", "Freq"], " moderate")
## [1] "According to the AMSTAR-2 scoring, 4 of these 31 reports were of high quality, 10 of low quality, and 17 of critically low quality.  moderate"
paste0("The interventions that covered the largest number of outcomes were ", paste(paste0(res_interventions$Intervention[1:5], " (n=", res_interventions$n[1:5], "/",length(unique(res_p$Outcome)), ")"), collapse=", "))
## [1] "The interventions that covered the largest number of outcomes were OXYT (n=11/19), PUFA (n=11/19), AAI (n=10/19), ACUP (n=9/19), NAC (n=9/19)"
# paste0("The interventions that covered the smallest number of outcomes were ", paste(paste0(res_interventions$Intervention[(nrowInter-8):nrowInter], " (n=", res_interventions$n[(nrowInter-8):nrowInter], "/",length(unique(res_p$Outcome)), ")"), collapse=", "))

res_outcomes = res_p %>%
  group_by(Outcome) %>%
  select(Intervention) %>%
  distinct() %>%
  summarise(n = n(),
            outcome = collapsunique(Intervention)) %>%
  arrange(desc(n))
## Adding missing grouping variables: `Outcome`
nrowOutcome = nrow(res_outcomes)
paste0("and the outcomes that were covered by the largest number of interventions all regarded the efficacy of the interventions on core ASD symptoms: ", paste(paste0(tolower(res_outcomes$Outcome[1:5]), " (n=", res_outcomes$n[1:5], "/",length(unique(res_m$Intervention)), ")"), collapse=", "))
## [1] "and the outcomes that were covered by the largest number of interventions all regarded the efficacy of the interventions on core ASD symptoms: social-communication (n=18/21), overall asd symptoms (n=16/21), restricted/repetitive behaviors (n=13/21), disruptive behaviors (n=12/21), adhd symptoms (n=9/21)"
res_adverse = res_p %>%
  filter(Outcome %in% safety) %>%
  group_by(Outcome, Intervention) %>%
  distinct() %>%
  summarise(n = n()) %>%
  group_by(Intervention) %>%
  summarise(n = n())
## `summarise()` has grouped output by 'Outcome'. You can override using the
## `.groups` argument.
paste0("In contrast, only ", nrow(res_adverse[res_adverse$n==3, ]),
       " interventions were tested for all our safety outcomes, namely acceptability, tolerability, and the presence of adverse effects (", paste(sort(unlist(res_adverse[res_adverse$n==3, "Intervention"])), collapse=", "), "). Some others were tested for one or two of these outcomes (", paste(sort(unlist(res_adverse[res_adverse$n!=3, "Intervention"])), collapse=", "), ").")
## [1] "In contrast, only 4 interventions were tested for all our safety outcomes, namely acceptability, tolerability, and the presence of adverse effects (MELAT, NAC, OXYT, PUFA). Some others were tested for one or two of these outcomes (L-CARNO, PROB, SULFO, V1a-RA, VIT-D)."
# sort(unique(lapply(res_outcomes[res_outcomes$Outcome %in% safety, "outcome"], function(x) unlist(strsplit(x, ", ")))[[1]]))
# res_outcomes[res_outcomes$Outcome %in% safety, ]
# sort(unique(lapply(res_outcomes[res_outcomes$Outcome %in% safety, "outcome"], function(x) unlist(strsplit(x, ", ")))[[1]]))



# paste0("The outcomes that were covered by the least number of interventions were ", paste(paste0(res_outcomes$Outcome[(nrowOutcome-5):nrowOutcome], " (n=", res_outcomes$n[(nrowOutcome-5):nrowOutcome], "/",length(unique(res_p$Intervention)), ")"), collapse=", "))

paste0("The median number of RCTs per meta-analysis was ",
      quantile(res_p$n_studies, 0.5), " (IQR = [", 
      quantile(res_p$n_studies, 0.25), ", ",
      quantile(res_p$n_studies, 0.75), "]), and the median number of participants per meta-analysis was ",
      quantile(res_p$total_n, 0.5), " (IQR = [", 
      quantile(res_p$total_n, 0.25), ", ",
      quantile(res_p$total_n, 0.75), "])")
## [1] "The median number of RCTs per meta-analysis was 4 (IQR = [2, 5]), and the median number of participants per meta-analysis was 154 (IQR = [92, 261.25])"
paste0("The median percentage of female per meta-analysis was ",
      round(quantile(res_p$perc_female, 0.5, na.rm=TRUE)), "% (IQR = [", 
      round(quantile(res_p$perc_female, 0.25, na.rm=TRUE)), ", ",
      round(quantile(res_p$perc_female, 0.75, na.rm=TRUE)), "])")
## [1] "The median percentage of female per meta-analysis was 14% (IQR = [10, 19])"
tab_Age = as.data.frame(table(res_p$Age))
paste0("A total of ",
       tab_Age[tab_Age$Var1=="< 6 yo", "Freq"], 
       " meta-analyses (",
       round(tab_Age[tab_Age$Var1=="< 6 yo", "Freq"] / nrow(res_p) * 100),
       "%) involved very young children (<6 years old), ", 
        tab_Age[tab_Age$Var1== "6-12 yo", "Freq"], 
       " meta-analyses (",
      round(tab_Age[tab_Age$Var1 == "6-12 yo", "Freq"] / nrow(res_p) * 100), "%) involved school-aged children (6-12 years old), and ",
        tab_Age[tab_Age$Var1== "13-19 yo", "Freq"], " meta-analyses (",
round(tab_Age[tab_Age$Var1 == "13-19 yo", "Freq"] / nrow(res_p) * 100), "%) involved adolescents (13-19 years old), and ",
        tab_Age[tab_Age$Var1== ">= 20 yo", "Freq"], " meta-analyses (",
round(tab_Age[tab_Age$Var1 == ">= 20 yo", "Freq"] / nrow(res_p) * 100), "%) involved adults (>= 20 years old)")
## [1] "A total of 26 meta-analyses (18%) involved very young children (<6 years old), 87 meta-analyses (59%) involved school-aged children (6-12 years old), and 15 meta-analyses (10%) involved adolescents (13-19 years old), and 20 meta-analyses (14%) involved adults (>= 20 years old)"
tab_IQ = as.data.frame(table(res_p$IQ))
paste0("A total of ",
       tab_IQ[tab_IQ$Var1=="Average (80-119)", "Freq"], " meta-analyses (",
       round(tab_IQ[tab_IQ$Var1=="Average (80-119)", "Freq"] / nrow(res_p) * 100), "%) focused on participants with overall cognitive functioning within the norm (IQ = [80-119]), ", 
        sum(tab_IQ[tab_IQ$Var1 %in% c("Limit (70-79)", "Low (< 70)"), "Freq"]), " meta-analyses (",
round(sum(tab_IQ[tab_IQ$Var1 %in% c("Limit (70-79)", "Low (< 70)"), "Freq"]) / nrow(res_p) * 100), "%) focused on participant with moderate or important cognitive difficulties (IQ < 80), and ",
        nrow(res_p) - sum(tab_IQ$Freq), " (",
       round(100-sum(tab_IQ$Freq) / nrow(res_p) * 100), "%) did not provide information on the cognitive functioning of the participants.")
## [1] "A total of 66 meta-analyses (45%) focused on participants with overall cognitive functioning within the norm (IQ = [80-119]), 23 meta-analyses (16%) focused on participant with moderate or important cognitive difficulties (IQ < 80), and 59 (40%) did not provide information on the cognitive functioning of the participants."
res_large = res_p %>% filter(eG > 0.5)
res_sig = res_p %>% filter(as.numeric(p_value) < 0.05)
res_effect = res_p %>% filter(eG > 0 & (eG > 0.5 | as.numeric(p_value) < 0.05))

paste0("The overall results indicated that ", 
       sum(res_p$eG > 0), " meta-analyses (", propC(res_p$eG > 0), 
       "%) had a pooled effect size in favor of the experimental group (SMD > 0, OR/RR > 1). More precisely, ", 
       nrow(res_large), " (",
       round(nrow(res_large) / nrow(res_p) * 100) , "%) had a magnitude that can be considered as moderate or large (SMD > 0.50), and ",
       nrow(res_sig), " (",
       round(nrow(res_sig) / nrow(res_p) * 100) , "%) had a statistically significant p-value (p < 0.05)")
## [1] "The overall results indicated that 113 meta-analyses (76%) had a pooled effect size in favor of the experimental group (SMD > 0, OR/RR > 1). More precisely, 26 (18%) had a magnitude that can be considered as moderate or large (SMD > 0.50), and 33 (22%) had a statistically significant p-value (p < 0.05)"
paste0("It should be acknowledged that the ",
       nrow(res_effect), " meta-analyses with moderate to large pooled effect sizes and/or statistically significant results are not without limitations, given that only ",
       sum(res_effect$rob > 75), " (", propC(res_effect$rob > 75), "%) of them had more than 75% of their participants coming from low-risk studies.",
       " Moreover, ",
       sum(res_effect$I2 >= 25), " (", propC(res_effect$I2 >= 25),
       "%) of these meta-analyses had a moderate or large inconsistency (I² >= 25%) and only ",
       sum(as.numeric(res_effect$PI_lo_g) > 0 | as.numeric(res_effect$PI_up_g) < 0, na.rm=TRUE), " (", 
       propC(as.numeric(res_effect$PI_lo_g) > 0 | as.numeric(res_effect$PI_up_g) < 0), "%) of the meta-analyses had a 95% prediction interval that excluded the null. Moreover, in ", propC(res_effect$n_agreement < 0.8), "% of the meta-analyses, more than 25% of the individual effect sizes were of the opposite sign to the pooled effect size.")
## [1] "It should be acknowledged that the 38 meta-analyses with moderate to large pooled effect sizes and/or statistically significant results are not without limitations, given that only 5 (13%) of them had more than 75% of their participants coming from low-risk studies. Moreover, 16 (42%) of these meta-analyses had a moderate or large inconsistency (I² >= 25%) and only 5 (13%) of the meta-analyses had a 95% prediction interval that excluded the null. Moreover, in NaN% of the meta-analyses, more than 25% of the individual effect sizes were of the opposite sign to the pooled effect size."
paste0(
  ifelse(nrow(res_p %>% filter(eG < 0  & as.numeric(p_value) < 0.05)), " there is a detrimental effect", "there is no detrimental effect" 
))
## [1] "there is no detrimental effect"

GRADING

nrow(res_p[res_p$GRADE == "High", ])
## [1] 0
dat_grade = as.data.frame(table(res_p$GRADE))
paste0("The assessment of the certainty of evidence allowed the emergence of a very clear pattern. Of the ", nrow(res_p), " unique PICOs retained in our main analysis, ", 
       dat_grade[dat_grade$Var1 == "High", "Freq"],
       " reached a 'High' GRADE rating, ", 
       dat_grade[dat_grade$Var1 == "Moderate", "Freq"],
       " reached a 'Moderate' rating, ",
       dat_grade[dat_grade$Var1 == "Low", "Freq"],
       " (", round(dat_grade[dat_grade$Var1 == "Low", "Freq"] / sum(dat_grade$Freq) * 100),"%) reached a 'Low' rating, and ",
       dat_grade[dat_grade$Var1 == "Very low", "Freq"],
       " (", round(dat_grade[dat_grade$Var1 == "Very low", "Freq"] / sum(dat_grade$Freq) * 100), "%) reached a 'Very low' rating. ")
## [1] "The assessment of the certainty of evidence allowed the emergence of a very clear pattern. Of the 148 unique PICOs retained in our main analysis, 0 reached a 'High' GRADE rating, 3 reached a 'Moderate' rating, 12 (8%) reached a 'Low' rating, and 133 (90%) reached a 'Very low' rating. "
low_mod_sig = dat_low_or_higher[dat_low_or_higher$eG > 0 & as.numeric(dat_low_or_higher$p_value) < 0.05, ]
paste0("More precisely, when considering only PICOs with a ‘Low’ or ‘Moderate’ rating, only ", nrow(low_mod_sig))
## [1] "More precisely, when considering only PICOs with a ‘Low’ or ‘Moderate’ rating, only 1"

Overlap

res_over <- res_m %>%
  group_by(PICO_amstar) %>%
  filter(n() > 1) %>%
  ungroup() %>%
  mutate(n_SRMA = length(unique(paper))) %>%
  group_by(PICO_amstar) %>%
  summarise(
    n = n(),
    n_SRMA = unique(n_SRMA),
    paper = collapsunique(paper),
    factor = collapsunique(Factor),
    min_es = min(eG),
    max_es = max(eG),
    dif_es = max(eG) - min(eG),
    max_grade = head(sort(GRADE), 1),
    min_grade = tail(sort(GRADE), 1),
    prop_sig = sum(as.numeric(p_value) < 0.05) / n(),
    avg_percentage_overlap = mean(map_dbl(combn(n(), 2, simplify = FALSE), function(idx) {
    interval_overlap_percentage(
      ci_lo_g[idx[1]], ci_up_g[idx[1]], 
      ci_lo_g[idx[2]], ci_up_g[idx[2]])
  })))

res_test = res_m %>%
 group_by(PICO_amstar) %>%
  filter(n() > 1) %>%
  mutate(PICO_hom = grepl("Homo", PICO_amstar_precise, fixed=TRUE)) %>%
  filter((PICO_hom & Age_precise == "Homogeneous") | !PICO_hom | IN_meta == 1) %>%
  group_by(PICO_amstar) %>%
  filter(n() > 1) %>%
  mutate(eG_meta_IN = eG[IN_meta == 1],
         GRADE_meta_IN = GRADE[IN_meta == 1],
         p_value_meta_IN = p_value[IN_meta == 1]) %>%
  ungroup() %>% rowwise() %>%
  filter(
    (abs(eG - eG_meta_IN) >= 0.50 & as.numeric(p_value) != as.numeric(p_value_meta_IN)) | 
    (GRADE == GRADE_meta_IN & as.numeric(p_value) != as.numeric(p_value_meta_IN)) |
    abs(as.numeric(GRADE) - as.numeric(GRADE_meta_IN)) >= 2 |
           IN_meta == 1) %>%
  # select(c(1:15, PICO_amstar, PICO_amstar_precise, IN_meta, amstar.x, Age, Age_precise, eG_meta_IN, GRADE, GRADE_meta_IN)) %>%
  arrange(PICO_amstar)
         
         
res_over_factors = res_test %>% 
  filter(abs(as.numeric(GRADE_meta_IN) - as.numeric(GRADE)) >= 2 | 
           (abs(eG_meta_IN - eG) >= 0.30 & sum(as.numeric(p_value_meta_IN) < 0.05,
                                               as.numeric(p_value) < 0.05) == 1))

res_over_disc = res_test %>%
  filter(Factor %in% res_over_factors$Factor |
         PICO_amstar %in% res_over_factors$PICO_amstar & IN_meta==1)
res_m$inter_out = paste0(res_m$Intervention, " - ", res_m$Outcome)
metaumbrella::forest(res_m %>% filter(PICO_amstar %in% res_over_disc$PICO_amstar),
       layout = "RevMan5",
       subgroup = "inter_out",
       subgroup.name = "",
       leftcols = c("paper", "Age", "Age_precise", "GRADE", 
                    "n_studies",  "rob_num", "active_controls", "I2", "effect.ci"),
       leftlabs = c("Paper", "Age", "Age_precise", "GRADE",
                    "n-studies", "Low RoB", "Active controls",  "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Key overlapping situations"
)

paste0("Among the ", nrow(res_over), " PICOs that were explored by at least 2 overlapping meta-analyses, the median overlap of the 95% CI of the pooled effect sizes was ",
       round(quantile(res_over$avg_percentage_overlap, 0.5)),
       "% (IQR=[", 
       round(quantile(res_over$avg_percentage_overlap, 0.25)),
       "%, ",
       round(quantile(res_over$avg_percentage_overlap, 0.75)),
       "%]). Only ", 
       sum(res_over$min_grade == res_over$max_grade), " (", 
       round(sum(res_over$min_grade == res_over$max_grade) / nrow(res_over) * 100), 
       "%) of PICOs consistently achieved a similar GRADE rating and only ", 
       sum(res_over$prop_sig %in% c(0,1)), " (", 
       round(sum(res_over$prop_sig %in% c(0,1)) / nrow(res_over) * 100), 
       "%) of PICOs consistently achieved a similar statistical significance status.") 
## [1] "Among the 53 PICOs that were explored by at least 2 overlapping meta-analyses, the median overlap of the 95% CI of the pooled effect sizes was 55% (IQR=[41%, 68%]). Only 36 (68%) of PICOs consistently achieved a similar GRADE rating and only 36 (68%) of PICOs consistently achieved a similar statistical significance status."
paste0("We found ", length(unique(res_over_disc$PICO_amstar)), " PICOs that presented a large discrepancy regarding the magnitude of the estimated effect (i.e., pooled effect sizes varying by more than SMD >= 0.30 with different statistical significance status), or regarding the confidence or GRADE status differing by two points (very low versus moderate; Fig X.). This information is displayed in details in Supplementary Materials SX, and we will illustrate this overlaping issue with the only PICO that generated ")
## [1] "We found 4 PICOs that presented a large discrepancy regarding the magnitude of the estimated effect (i.e., pooled effect sizes varying by more than SMD >= 0.30 with different statistical significance status), or regarding the confidence or GRADE status differing by two points (very low versus moderate; Fig X.). This information is displayed in details in Supplementary Materials SX, and we will illustrate this overlaping issue with the only PICO that generated "

OXYT restrict

dat_oxyt_rest = dcct %>% 
   filter(intervention_general == "NAC" &
          first_author_meta %in% c("Iffland", "Siafis (child)") &
          outcome_general == "Disruptive behaviors") %>%
  arrange(author, year) %>%
  mutate(TE = ifelse(!is.na(reverse_es), -as.numeric(as.character(value)), value),
         TE_lo = ifelse(!is.na(reverse_es), -as.numeric(as.character(ci_up)), ci_lo),
         TE_up = ifelse(!is.na(reverse_es), -as.numeric(as.character(ci_lo)), ci_up), 
         ID = paste0(author, " (", year, ")"))

met_oxyt_res = meta::metagen(TE = as.numeric(TE),
   lower = as.numeric(TE_lo),
   upper = as.numeric(TE_up),
   subgroup = ID,
   subgroup.name = "Study",
   studlab = meta_review,
   data = dat_oxyt_rest)


m2 <- meta::rob(RoB_SeqGen, RoB_Conceal, RoB_BlindPers, RoB_BlindOut, RoB_Attrition, RoB_Reporting, overall = rob, data = met_oxyt_res, tool = "rob1")

meta::forest(m2,
       layout = "RevMan5", common = FALSE, random = FALSE, 
       leftcols = c("studlab", "effect.ci"),
       leftlabs = c("Study","SMD [95% CI]"),
       col.square = "gray", 
  col.study = "black", col.square.lines = "black",
  overall = FALSE, hetstat = FALSE)

#        subgroup = "author",
#        subgroup.name = ""
# metaumbrella::forest(focus_oxyt_restrict_adult,
#        layout = "RevMan5",
#        subgroup = "inter_out",
#        subgroup.name = "",
#        leftcols = c("paper", "Age", "Age_precise", "GRADE", 
#                     "n_studies",  "rob", "active_controls", "I2", "effect.ci"),
#        leftlabs = c("Paper", "Age", "Age_precise", "GRADE",
#                     "n-studies", "Low RoB", "Active controls",  "I²", "eSMD + 95% CI"),
#        # sortvar = res_p$Age,
#        smlab = "Key overlapping situations"
# )
res_over_sig <- res_m %>%
  group_by(PICO_amstar) %>%
  filter(n() > 1) %>%
  mutate(eG_meta_IN = eG[IN_meta == 1],
         GRADE_meta_IN = GRADE[IN_meta == 1],
         p_value_meta_IN = p_value[IN_meta == 1]) %>%
  filter(as.numeric(p_value_meta_IN) < 0.05) %>%
  summarise(
    n = n(),
    paper = collapsunique(paper),
    factor = collapsunique(Factor),
    min_es = min(eG),
    max_es = max(eG),
    dif_es = max(eG) - min(eG),
    max_grade = head(sort(GRADE), 1),
    min_grade = tail(sort(GRADE), 1),
    prop_sig = sum(as.numeric(p_value) < 0.05) / n(),
    avg_percentage_overlap = mean(map_dbl(combn(n(), 2, simplify = FALSE), function(idx) {
    interval_overlap_percentage(
      ci_lo_g[idx[1]], ci_up_g[idx[1]], 
      ci_lo_g[idx[2]], ci_up_g[idx[2]])
  })))


nrow(res_over_sig %>% filter(prop_sig < 1)) / nrow(res_over_sig)
## [1] 0.6
res_discrepant = res_m %>% filter(PICO_amstar %in% res_over_disc$PICO_amstar)
metaumbrella::forest(res_discrepant,
       layout = "RevMan5",
       subgroup = "PICO_amstar",
       subgroup.name = "",
       leftcols = c("Intervention", "Outcome", "Age_precise", "paper", 
                    "GRADE", 
                    "n_studies",  "rob_num",  
                    "I2", "effect.ci"),
       leftlabs = c("Intervention", "Outcome", "Homogeneous Age", "Paper", 
                    "GRADE",
                    "n-studies", "Low RoB",   
                    "I²", "eSMD + 95% CI"),
       # sortvar = res_p$Age,
       smlab = "Discrepant PICOs"
)

propC((res_over$min_grade == "Very low" & res_over$max_grade == "Moderate") |
       res_over$dif_es>= 0.50 & !res_over$prop_sig %in% c(0, 1))
## [1] 15
paste0("The median overlap per meta-analysis was ",
      quantile(res_over$avg_percentage_overlap, 0.5), " (IQR = [", 
      quantile(res_over$avg_percentage_overlap, 0.25), ", ",
      quantile(res_over$avg_percentage_overlap, 0.75), "])")
## [1] "The median overlap per meta-analysis was 55.394641564084 (IQR = [41.0098522167488, 68.2401532502682])"
average_consistency <- dat %>%
  group_by(intervention) %>%
  summarise(avg_percentage_overlap = mean(map_dbl(combn(n(), 2, simplify = FALSE), function(idx) {
    interval_overlap_percentage(
      lower_bound[idx[1]], upper_bound[idx[1]], 
                                
      lower_bound[idx[2]], upper_bound[idx[2]])
  })))
res_m_cred = metaumbrella::add.evidence(res,
                           criteria = "Personalized",
                           class_I = c(total_n = 500, p_value = 1e-3, I2 = 50,
                                       egger_p = .05, pi = "notnull",
                                       rob = 75),
                           class_II = c(total_n = 350, p_value = 1e-3,
                                        largest_CI = "notnull",
                                        rob = 50),
                           class_III = c(total_n = 200, p_value = 1e-3),
                           class_IV = c(p_value = 5e-2), verbose = FALSE)